Pearson Chi-squared test with approximated methods for residual degrees of freedom in GLMMs
Explanations
package: glmmrBase https://cran.r-project.org/web/packages/glmmrBase/index.html
website (out of date): https://samuel-watson.github.io/glmmr-web/
preprint about the package (2024): https://arxiv.org/abs/2303.12657
mantainer: Samuel Watson - Prof stats Uni Birmingham https://www.birmingham.ac.uk/staff/profiles/applied-health/watson-sam
To create the Pearson Chi-squared test with the approximate residual degrees of freedom, I first:
- converted the lme4 function to the glmmrBase function
- and used the method small-sample_correction() to have the dfs
correction:
- Kenward-Roger
- Kenward-Roger2
- Satterthwaithe
- Kenward-Roger
- I subtracted the sample size from this df to have the residual
df.
- I used two.sides test.
From Watson 2024: “The degrees of freedom correction is also given in Kenward and Roger (1997), which is the same degrees of freedom correction originally proposed by Satterthwaite (1946). The Kenward-Roger correction, though, can over-estimate the standard errors in very small samples in some cases. Kenward and Roger (2009) proposed an improved correction for covariance functions not linear in parameters (such as the autoregressive or exponential functions). The improved correction adds an additional adjustment factor. We have implemented these corrections in the package, which can returned directly from a Model object using the small_sample_correction() member function”
The function:
approximateDFpearson <- function(model, data, type=c("naive", "KR", "KR2", "sat"),
alternative=c("two.sided", "greater", "less")){
if(class(model)[1] != "glmerMod") stop("model is not a glmerMod from lme4")
if(model@resp$family$family != "poisson") stop("model is not a Poisson")
if(type == "naive"){
rdf <- df.residual(model)
} else{
f1 <- glmmrBase::lme4_to_glmmr(formula(model), colnames(data))
mod <- glmmrBase::Model$new(f1, data=data, family=poisson())
rdf <- dim(data)[1]- mod$small_sample_correction(type)$dof[1]
}
rp <- residuals(model, "pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
if(alternative == "greater") pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
else if (alternative == "less") pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE)
else if (alternative == "two.sided") pval <- min(min(pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE), pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)) * 2,1)
out = list()
out$statistic = prat
names(out$statistic) = "dispersion"
out$parameter = rdf
names(out$parameter) = paste0("df-",type)
out$method = "Parametric dispersion test via mean Pearson-chisq statistic"
out$alternative = alternative
out$p.value = pval
class(out) = "htest"
return(out)
}
Loading simulation outputs
Simulations done in 8_glmm_approxiDFpearson.R:
Sim setting:
sampleSize <- c(200, 500, 1000)
intercept <- c(-1.5,0,1.5)
overdispersion = seq(0,1,0.1)
ngroups <- c(10, 50, 100)
nRep = 1000
slope = 1
randomEffectVariance = 1load(here("data", "8_approximateDFpearson.Rdata")) # simulated data
(tspent <- map_dbl(out.pois, "time")) # minutes## 10_200_-1.5 10_200_0 10_200_1.5 10_500_-1.5 10_500_0
## 1.321605 1.219799 1.216456 2.353973 2.326321
## 10_500_1.5 10_1000_-1.5 10_1000_0 10_1000_1.5 50_200_-1.5
## 2.370980 8.441990 8.488686 8.568489 1.519092
## 50_200_0 50_200_1.5 50_500_-1.5 50_500_0 50_500_1.5
## 1.471159 1.405490 3.685771 6.909883 3.672749
## 50_1000_-1.5 50_1000_0 50_1000_1.5 100_200_-1.5 100_200_0
## 16.131561 16.184664 16.255666 1.805094 1.823181
## 100_200_1.5 100_500_-1.5 100_500_0 100_500_1.5 100_1000_-1.5
## 1.689025 5.617923 5.463803 5.451093 28.755045
## 100_1000_0 100_1000_1.5
## 31.243332 28.750171
#simulations
simuls.pois <- map_dfr(out.pois, "simulations", .id="model") %>%
separate(model, c("ngroups", "sampleSize", "intercept"), sep="_") %>%
rename(overdispersion = controlValues,
PearNaive.rdf = PearNaive.rdf.df.naive,
PearSat.rdf = PearSat.rdf.df.sat,
PearKR.rdf = PearKR.rdf.df.KR,
PearKR2.rdf = PearKR2.rdf.df.KR2)There were many negative residual dfs, it generates NaN to the test:
It was around 6% of the simulations.
Almost everything came from the 100 groups in 200 sample sizes - 2 observations per group:
##
## 200 500
## 100 15567 392
## 50 1572 0
The percentage of failure was around the same among the corrections
## [1] 0.02415152
## [1] 0.02332997
## [1] 0.02317508
Subsetting the simulation output only for the positive rdfs:
Comparing corrections of residual df
dof <- simuls.pois2 %>% dplyr::select(contains("rdf"), replicate, ngroups,
overdispersion, intercept, sampleSize) %>%
pivot_longer(1:4, names_to = "test", values_to = "res.df") %>%
group_by(sampleSize, ngroups, intercept,overdispersion, test)
dof$intercept <- fct_relevel(dof$intercept, "-1.5", "0", "1.5")
dof$ngroups <- fct_relevel(dof$ngroups, "10", "50", "100")
dof$sampleSize <- as.factor(as.numeric(dof$sampleSize))Very similar rdfs. Satterthwaite had a bit smaller rdfs for small
sample size.
Type I error
p.pois <- simuls.pois2 %>% dplyr::select(ends_with(".p"), replicate,
ngroups,
overdispersion, intercept, sampleSize) %>%
pivot_longer(1:4, names_to = "test", values_to = "p.val") %>%
group_by(sampleSize, ngroups, intercept, overdispersion, test) %>%
summarise(p.sig = sum(p.val<0.05,na.rm=T),
nsim = length(p.val[!is.na(p.val)]))## `summarise()` has grouped output by 'sampleSize', 'ngroups', 'intercept',
## 'overdispersion'. You can override using the `.groups` argument.
p.pois$prop.sig <- p.pois$p.sig/p.pois$nsim
p.pois$intercept <- fct_relevel(p.pois$intercept, "-1.5", "0", "1.5")
p.pois$ngroups <- fct_relevel(p.pois$ngroups, "10", "50", "100")
p.pois$sampleSize <- as.factor(as.numeric(p.pois$sampleSize))Power
- Power for the 3 correction methods is basically the same.
- Corrected power still very bad for small intercepts (-1.5 and 0),
figure for each N groups separated
Dispersion statistics
Median of disp stat
d.pois <- simuls.pois2 %>% dplyr::select(ends_with("dispersion"), replicate, ngroups,
overdispersion, intercept, sampleSize) %>%
pivot_longer(1:4, names_to = "test", values_to = "dispersion") %>%
group_by(sampleSize, ngroups, intercept,overdispersion, test) %>%
summarise(median.stat = median(dispersion, na.rm=T))## `summarise()` has grouped output by 'sampleSize', 'ngroups', 'intercept',
## 'overdispersion'. You can override using the `.groups` argument.
d.pois$intercept <- fct_relevel(d.pois$intercept, "-1.5", "0", "1.5")
d.pois$ngroups <- fct_relevel(d.pois$ngroups, "10", "50", "100")
d.pois$sampleSize <- as.factor(as.numeric(d.pois$sampleSize))Simplified figure with just one correction (KR2) to compare with the
Naïve
Compare power with the other tests
Compare pearson corrections with Dharma-conditional and pearson par boot.
Using only ONE correction, 100 groups
load(here("data", "5_glmmPois_power_50.Rdata")) # simulated data
out.pois50 <- out.pois
load(here("data", "5_glmmPois_power_100.Rdata")) # simulated data
out.pois100 <- out.pois
out.pois <- flatten(list(out.pois50, out.pois100))
simuls.pois <- map_dfr(out.pois, "simulations", .id="model") %>%
separate(model, c("ngroups", "sampleSize", "intercept"), sep="_") %>%
rename("overdispersion" = "controlValues")
# power
pw <- simuls.pois %>% dplyr::select(dhaCO.p.val,refCO.p.val, replicate,
ngroups,
overdispersion, intercept, sampleSize) %>%
filter(intercept %in% c( "-1.5", "0", "1.5"),
sampleSize %in% c("200", "500", "1000")) %>%
pivot_longer(1:2, names_to = "test", values_to = "p.val") %>%
group_by(sampleSize, ngroups, intercept, overdispersion, test) %>%
summarise(p.sig = sum(p.val<0.05,na.rm=T),
nsim = length(p.val[!is.na(p.val)]))## `summarise()` has grouped output by 'sampleSize', 'ngroups', 'intercept',
## 'overdispersion'. You can override using the `.groups` argument.
pw$prop.sig <- pw$p.sig/pw$nsim
pw$intercept <- fct_relevel(pw$intercept, "-1.5", "0", "1.5")
pw$sampleSize <- as.factor(as.numeric(pw$sampleSize))comb <- p.pois %>% filter(ngroups %in% c("50", "100"), test == "PearKR2.p") %>%
bind_rows(pw)
ggplot(comb, aes(x=overdispersion, y=prop.sig, col=test, linetype=ngroups, shape=ngroups))+
geom_point(alpha=0.7) + geom_line(alpha=0.7) +
facet_grid(sampleSize~intercept) +
geom_hline(yintercept = c(0.05,0.5), linetype="dotted") +
ylim(0,1)+
theme(panel.background = element_rect(color="black"),
legend.position = "bottom") +
guides(color=guide_legend(nrow=3, byrow=TRUE),
linetype=guide_legend(nrow=2, byrow=TRUE),
shape=guide_legend(nrow=2, byrow=TRUE)) +
xlab("Overdispersion") + ylab("Power")+
scale_color_discrete(labels = c("Sim-based conditional",
"Pearson Chi-sq KR2",
"Pearson Param. Boot."))